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ABSTRACT 

The recurrent outbursts that characterise low-mass binary systems reflect thermal state 
changes in their associated accretion discs. The observed outbursts are connected to the strong 
variation in disc opacity as hydrogen ionises near 5000 K. This physics leads to accretion disc 
models that exhibit bistability and thermal limit cycles, whereby the disc jumps between a 
family of cool and low accreting states and a family of hot and efficiently accreting states. 
Previous models have parametrised the disc turbulence via an alpha (or 'eddy') viscosity. In 
this paper we treat the turbulence more realistically via a suite of numerical simulations of 
the magnetorotational instability (MRI) in local geometry. Radiative cooling is included via 
a simple but physically motivated prescription. We show the existence of bistable equilibria 
and thus the prospect of thermal limit cycles, and in so doing demonstrate that MRl-induced 
turbulence is compatible with the classical theory. Our simulations also show that the tur- 
bulent stress and pressure perturbations are only weakly dependent on each other on orbital 
times; as a consequence, thermal instability connected to variations in turbulent heating (as 
opposed to radiative cooling) is unlikely to operate, in agreement with previous numerical 
results. Our work presents a first step towards unifying simulations of full MHD turbulence 
with the correct thermal and radiative physics of the outbursting discs associated with dwarf 
novae, low-mass X-ray binaries, and possibly young stellar objects. 
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1 INTRODUCTION 

Recurrent outbursts in accreting systems are commonly attributed 
to global instabilities in their associated accretion discs. In particu- 
lar, it is believed that thermal instability driven by strong variations 
in the disc's cooling rate causes the observed state transitions in 
dwarf novae (DNe) and low-mass X-ray binaries (LMXBs) (La- 
sota 2001), while the rich array of accretion variability associated 
with quasars could be excited by thermal instabiUty driven by vari- 
ations in the turbulent heating rate (Shakura and Sunyaev 1976, 
Abramowicz et al. 1988), by thermal-viscous instability (Light- 
man & Eardley 1974), or assorted dynamical instabilities (e.g. Pa- 
paloizou & Pringle 1984, Kato 2003, Ferreira & Ogilvie 2009). 
On the other hand, FU Ori outbursts, characteristic of protostellar 
discs, probably result from the interplay of thermal, gravitational, 
and magnetorotational instability (MRI) across the intermittently 
inert dead zone (Gammie 1996, Balbus & Hawley 1998, Armitage 
et al. 2002, Zhu et al. 2009, 2010). Exceptions to this class of model 
include classical novae whose outbursts can be traced to thermonu- 
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cleur fusion of accreted material on the white dwarf surface (Gal- 
lagher & Starrfield 1978, Shara 1989, Starrfield et al. 2000). 

The most developed, and possibly most successful, model of 
accretion disc outbursts describes the recurrent eruptions that char- 
acterise the light curves of DNe and LMXBs. In this model, the 
ionisation of hydrogen, occurring at temperatures of around 5000 
K, induces a significant opacity change in the disc orbiting the pri- 
mary star (Faulkner et al. 1983, hereafter FLP83), which then leads 
to thermal instability and hysteresis in the disc gas (Pringle 1981). 
The system exhibits a characteristic 'S-curve' in the phase plane 
of surface density E and central temperature Tc (cf. Fig. 1), and 
as a consequence the gas falls into one of two stable states: an 
optically thick hot state, characterized by strong accretion, or an 
optically thin cooler state, in which accretion is less efficient (e.g. 
Meyer & Meyer-Hofmeister 1981, FLP83). Outbursts can then be 
modelled via a 'limit cycle', whereby the disc jumps from the low 
accreting state to the high accreting state and then back again as 
mass builds up and is then evacuated. This mechanism, based on 
local bistable equilibria, is the foundation for a variety of advanced 
models, which have enjoyed significant successes in reproducing 
the observed behaviour of accretion discs in binary systems, even 
if interesting discrepancies persist (Lasota 2001). 
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The classical model of DNe and LMXBs treats the disc as 
laminar and assumes that the disc turbulence can be modelled with 
the a prescription (Shakura & Sunyaev 1973), whereby the action 
of the turbulent stresses is characterized as a diffusive process with 
an associated 'eddy' viscosity (Balbus & Papaloizou 1999). This 
has been sufficient to sketch out the qualitative features of puta- 
tive outbursts, but such a crude description has imposed limitations 
on the formalism that are now impeding further progress (Lasota 
2012). On the other hand, fully consistent magnetohydrodynamic 
(MHD) simulations of disc turbulence generated by the MRI have 
been performed for almost 20 years (Hawley et al. 1995, Stone et 
al. 1996, Hawley 2000, etc), though typically with simplified ther- 
modynamics (e.g. isothermaUty). It is the task of this paper to begin 
the process of uniting the classic thermal instability models of DNe 
and LMXBs with full MHD simulations of the MRI, and thus con- 
sistently account for both the turbulence and radiative cooling. The 
first step we take is limited to local models, as even though discs 
undergo global outbursts, the ability of local annuli to exhibit hys- 
teresis behaviour is key. Local studies will permit us to assess if 
and how the classic model of thermal instability can operate in the 
presence of realistic turbulent heating. At the same time, they pro- 
vide an excellent test of the MRI itself; if the MRI is to remain the 
chief candidate driving disc accretion it must fulfil its obUgations 
to classical disc theory. 

We undertake unstratified shearing box simulations of the 
MRI that include Ohmic and viscous heating and a radiative cool- 
ing prescription that is able to mimic the transition between the 
optically thin and thick states (FLP83). These simulations clearly 
reproduce S-curves in the E and mean temperature plane, and these 
constitute the main result of our paper. We can thus animate local 
thermal limit cycles via a sequence of local box runs. In particular, 
the simulated turbulent heating is found to be 'well-behaved' and 
not so intermittent as to prematurely disrupt the thermal limit cycles 
required by the classical theory. In addition, simulations with dif- 
fering net toroidal and vertical fiuxes produce S-curves that exhibit 
variable mean values of a, suggesting that when global simulations 
are considered, and net-flux is no longer conserved locally, varia- 
tions in effective a between the upper and lower branches of the 
S -curve may be produced, as required by the classical model (e.g. 
Smak 1984a). 

Finally, the simulated short-term behaviours of the average 
viscous stress and the disc pressure reveal only a weak func- 
tional dependence, as remarked upon by previous authors (Hirose 
et al. 2009). This poor correlation means that proposed thermal in- 
stabiUties driven by a direct heating response to imposed pressure 
perturbations are imUkely to function in radiation-pressure domi- 
nated accretion discs (e.g. Shakura & Sunyaev 1976, Abramowicz 
et al. 1988). This is in marked contrast to the thermal instabilities 
implicated in DNe and LMXBs, which are driven by the cooling 
response to imposed temperature perturbations, mediated via opac- 
ity variations. Note, however, that over long times we find that the 
stress and pressure are correlated. 

The plan of the paper is as follows. In the following section 
we discuss the thermal instability model for DNc and LMXBs in 
more detail. In Section 3 we present the governing equations while 
outlining the radiative cooling prescription the simulations adopt. 
In section 4 wc give the numerical details of the simulations. Our 
results are presented in Section 5, and we conclude in Section 6. 



2 BACKGROUND 

We consider close binary systems that consist of one low-mass 
lobe-filling star transferring material to a degenerate companion, 
either a white dwarf or a black-hole/neutron star. The characteristic 
feature of these systems is their recurrent eruptive activity in the 
optical or X-ray spectnmi, respectively. DNe undergo outbursts of 
some 2-5 magnitude which last 2-20 days, with recurrence intervals 
ranging from 10 days to years (see e.g. Warner 1995). In addition, 
certain sources exhibit more complex behaviour, such as 'super- 
outbursts' and 'standstills' (whereby the cycle is interrupted for an 
indefinite period). LMXBs are poorly constrained relative to DNe 
because they emit almost exclusively in X-rays, and so do not enjoy 
the same observational coverage. Even so, it has been established 
that their outbursts exhibit luminosity enhancements of several or- 
ders of magnitude, while their spectra progress through a series of 
canonical X-ray states (McClintock and Remillard 2006). 

It is generally accepted that outbursts in these systems take 
place in the accretion disc that orbits around the primary star and 
are triggered once sufficient mass has built up in the disc. During 
an outburst the disc jumps from a cool low-accreting state to a hot 
high-accreting state that deposits this surfeit of material onto the 
central object. After this mass has been evacuated the disc returns 
to its low state and the cycle repeats. The basic outline of the model 
was first proposed by Smak (1971) and Osaki (1974); but the physi- 
cal basis for instability in the disc was not understood till much later 
(FLP83, see also Hoshi 1979). Researchers now recognise that the 
cycles are the consequence of a thermal instability related to the 
variation of opacity with temperature. 

Above roughly 5000K the coUisional ionisation of hydrogen 
commences, and this liberates an increasing number of free elec- 
trons that enhance the gas opacity k. The dependence of opacity on 
temperature, for a range of densities, is conveniently illustrated in 
Fig. 9 of Bell & Lin (1994). The figure shows that, after ionisation 
of hydrogen is complete T > lO'^K, we have roughly k oc T^^'^ , 
and thus k decreases with temperature. On the other hand, cold 
molecular gas, with T < 2 x 10^ K, exhibits a k that also decreases 
with temperature, primarily on account of dust grain destruction. 
But the gas opacity in the intermediate temperature regime be- 
tween these two limits increases rapidly with temperature, because 
of the fiood of free electrons. Here k oc T^", approximately. This 
dramatic tendency for heat to be better trapped as temperature in- 
creases leads immediately to thermal instabiUty in this regime. 

The instability criterion can be formulated quantitatively via 
the following argument. Let us denote by cd the local turbulent 
dissipation rate per unit area. Then local thermal equilibrium is de- 
termined by balancing this injection rate against the disc's radiative 
losses: 

ED = 2aT^, (1) 

where Te is the effective temperature of the disc's photosphere and 
a is Stefan's constant. If we now denote by Tc the midplane tem- 
perature of the disc, a local thermal stability criterion for this state 
is simply: 

(Lasota 2001). This expresses the requirement that when the mid- 
plane temperature Tc is increased, the cooling response via radia- 
tion must outstrip the increase in the heating. If this condition is not 
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Figure 1. Representative plot of an S-curve defining local steady state ther- 
mal equilibria with limit cycle arrows. This is constructed from the laminar 
a disc modelling of FLP83. For reference the parameters were E = 3.5, 
A = 5, /i = 1, and ags = 0.0275 (see Section 3 and the Appendix). We 
remark that in this figure, central temperature is plotted as ordinate. How- 
ever, if the mass flow rate M is used instead, a curve of the same qualitative 
form results (FLP83). 



agates through the disc (Papaloizou et al. 1983, Meyer & Meyer- 
Hofmeister 1984, Smak 1984a). Provided conditions are favourable 
and the front propagates inwards, converting the entire disc into an 
upper state, there will be a global outburst that results in a sig- 
nificant mass deposit onto the primary. Consequently, the surface 
density everywhere decreases until the disc undergoes a downward 
transition to the lower branch. 

Our discussion of the thermal equilibria and limit cycles need 
only consider local models of an accretion disc. However, more ad- 
vanced DN and LMXB models must incorporate additional global 
physics in order to capture the heating and cooling fronts that me- 
diate branch changes, and ultimately to match specific observations 
in detail (e.g. Papaloizou et al. 1983, Smak 1984b, Ichikawa & Os- 
aki 1992, Menou et al. 1999). Even so, all such global treatments 
are founded on the hysteresis behaviour of the local model (sin- 
gle annulus) outlined above. In this paper we concentrate on this 
fundamental engine of instability. Our aim is to show that it can 
function in the presence of turbulence generated by the MRI. This 
is intended as a first step in moving disc outburst modelling away 
from the heuristic a prescription. Future work can then include the 
detailed global physics. 



met, then the midplane temperature increase will be reinforced and 
there will be a thermal runaway. 

At high temperatures, well above the ionisation threshold, we 
find that the disc radiates efficiently and that dhiT* / dlaTc ~ 7. 
The heating rate e_D , on the other hand, increases much more slowly 
(oc Tc for a laminar a model) and so thermal stability ensues. How- 
ever, in the transition regime, where 4 x 10"^ < T < 10* K, the 
photospheric temperature varies very weakly with Tc because 
the opacity is increasing rapidly. The cooling rate barely responds 
to an increase in central temperature Tc, which leads to the trap- 
ping of excess heat and thermal instability. At the lowest temper- 
atures, T < 4000K, stability is recovered because the disc enters 
an optically thin regime in which the rate of surface radiation in- 
creases rapidly with Tc. In summary, the disc is bistable and will 
tend to move to one of two thermally stable steady regimes: (a) 
hot and ionised (hence efficiently accreting), with T > 10*K; or 
(b) cold and predominantly neutral (less efficiently accreting) with 
T < 4000K. In addition, there exists an unstable intermediate equi- 
librium state between these two limits. 

Because for a given surface density E there can be up to three 
Tc associated with a thermal equilibrium, the family of equilibrium 
solutions sketches out a characteristic S-curve in either the (E, Tc) 
plane or, alternatively, the (E, M) plane, where M is mass accre- 
tion rate. Note that M can be directly related to Tc and E (e.g. see 
FLP83 and section l42l below) . In Fig. I, we plot a representative 
family of such thermal equilibria in the (S, Tc) plane computed via 
the techniques of FLP83. There is a range of S for which the sys- 
tem supports three states, two stable (indicated with a blue colour), 
one unstable (indicated with a red colour). 

In Fig. I the trajectory of the outbursting limit cycle is repre- 
sented by the black arrows. Consider gas physically located at the 
outermost disc radii and in a thermodynamic state associated with 
the lower branch. At this radius, mass steadily accumulates because 
accretion onto the primary star cannot match the mass transfer into 
the disc from the secondary. Consequently, the surface density in- 
creases and we travel up the lower branch until it ends, at which 
point the gas heats up dramatically. Once it settles on the upper 
branch accretion is much more efficient and a transition wave prop- 



3 GOVERNING EQUATIONS 

We adopt the local shearing box model (Goldreich & Lynden-Bell, 
1965) and a Cartesian coordinate system (x, y, z) = {xi, X2, xs) 
with origin at the centre of the box. These coordinates represent the 
radial, azimuthal, and vertical directions respectively. Vertical strat- 
ification is neglected. The basic equations express the conservation 
of mass, momentum, and energy and include the induction equa- 
tion for the magnetic field. They incorporate a constant kinematic 
viscosity i/ and magnetic diffusivity 77. The continuity equation, 
equation of motion, and induction equations are given by 



| + V.(P.) = 0, 



{VxB)xB 

4tt 



VP + v-n, 



OB 

— = V X (1; X B - 77 V X B) 

and the first law of thermodynamics is written in the form 



P 



De P Dp 



n : Vi; + 



77|VxB| 



Dt p Dt . . - . 

Here the convective derivative is defined through 

D_ ^ d_ 

Dt " dt 



A. 



- V ■ V, 



(3) 

(4) 
(5) 

(6) 

(7) 



is the Keplerian orbital frequency evaluated at the centre of the 
box, V is the velocity, B is the magnetic field and e is the internal 
energy per unit mass. The tidal potential is "l? = —30,^x^/2 and 
A is the (radiative) cooling rate per unit volume. The right hand 
side of Equation ^ includes contributions from viscous dissipa- 
tion, Ohmic heating and radiative losses. The components of the 
viscous stress tensor 11 are given by 



, ' dvi dvk 2 
H,fc = pu [ h -difc V-v 

, (^Xk OXi o 



(8) 
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In addition, the magnetic field must satisfy V • i? = 0. Finally, we 
adopt an ideal gas equation of state 

P^'^P^pcl, (9) 

where — TZT / ^ is the isothermal sound speed, with T the tem- 
perature, TZ the gas constant, and fj, the mean molecular weight. 

3.1 Radiative cooling model 

In order to proceed, we need to realistically account for the radia- 
tive cooling of the gas, via the term A. The shearing box by defi- 
nition is located far from the disc's upper and lower surfaces with 
periodic boundary conditions applied in the vertical direction. Nev- 
ertheless, we can approximate the box's radiative losses in a physi- 
cally meaningful way via a prescription outlined in FLP83. 

Our model assumes that the cooling A is a function only of 
time and horizontal position [x, y). It is taken to be 



A = 2aT^/Ho, 



(10) 



where Te is interpreted as the effective temperature at the upper and 
lower vertical boundary of the disc. Here Hq is a reference scale 
height. The main task is to derive how the surface temperature 
varies in response to changes in the midplane temperature Tc- 

We envisage that Te takes a different form in three distinct 
physical disc regime^ 

(i) The first regime corresponds to the hot optically thick con- 
ditions associated with the disc's high accreting state. In such a 
state, the disc's photospheric temperature Te is marginally above 
the value appropriate for the ionisation threshold for hydrogen 
(~ 5000 K) and so the entire disc is optically thick. The classi- 
cal relationship — (4/3)(r^/rc), where Tc is midplane optical 
depth, can be used to relate Te and Te . As most of the optical depth 
comes from the regions close to the midplane, Te can be evaluated 
using the midplane opacity. 

(ii) The second regime corresponds to an intermediate 'hybrid' 
state in which the midplane is hot and partially ionised but the sur- 
face layers have become sufficiently cool for hydrogen ionisation 
to fall off sharply. As a consequence, the opacity drops near the 
photosphere with drastic consequences for the disc's vertical tem- 
perature structure, as in cool stars (Hayashi & Hoshi 1961). The 
classical dependence of Te on Te breaks down, with the coimection 
between the two becoming quite weak. Instead, we must determine 
Te from the outer boundary condition Te = 1, where Te is photo- 
spheric optical depth. 

(iii) The third regime corresponds to the cool optically thin 
regime in which the entire disc, including the midplane, lies 
marginally below the ionisation threshold. In this regime, consid- 
eration of direct cooling then yields a simple approximate relation- 
ship between Te and Te- 

In the Appendix we consider each regime in turn and obtain 
three functional forms relating Te to Te and E. Our final expression 
for Te is an interpolation formula, Eq. ( IA13t . connecting these 



^ FLP83 originally introduced 4 regimes, with an extra regime correspond- 
ing to the bottom right 'comer' of the S-curve where the gas may be opti- 
cally thick. It is in fact unnecessary to treat this regime separately, as it is 
covered by regimes 2 and 3. 



different regimes. In sunmiary, the prescription for finding Te is 
given by 



Te = < 



Nl/4 



(4/3r, , 
(^lO^** Epe^'^/T>i 



1/10 



' T 

^ e 



in Regime 1 , 
in Regime 2, 
in Regime 3, 



(11) 



where E and A are two dimensionless constants, and the midplane 
optical depth is Te — KcE with Kc being the opacity evaluated 
at the midplane density pe and temperature Te. All dimensional 
quantities are in cgs units. As in FLP83 the opacity in regime 1 is 
taken to be 



K = 1.5 X lO^Vr^^'^ 



and in regime 3 it is taken to be 

1^-30 l/SrrilO 



(12) 



(13) 



(14) 



In addition, 

E _En/^i\i/2 

Finally, we specify a constant p for all regimes. This of course ne- 
glects partial ionization, but should not interfere especially in our 
qualitative results. 

In FLP83, thermal equilibrium solutions are computed by 
matching the cooling rate, computed according to the above proce- 
dure, with the turbulent alpha heating, for which the volume heating 
rate is 



ED = 3assp7^TO/(2^i), 



(15) 



where the Shakura-Sunyaev ass parameter is a dimensionless con- 
stant less than one (cf. Eq. ll21t). These solutions may be viewed as 
relations between E and Te, though they are more commonly re- 
expressed as relations between E and Te . Some of these S-curves 
are illustrated in Figs[T]and 4-6. 



4 NUMERICAL SET-UP 

We solve equations (O-® with a parallel version of ZEUS (Stone 
& Norman 1992a, 1992b). This uses a finite difference scheme to 
obtain the spatial derivatives, is first order explicit in time, and em- 
ploys constrained transport to ensure B remains solenoidal. Our 
version of ZEUS has been altered according to improvements out- 
lined by Lesaffre & Balbus (2007), Silvers (2008), and Lesaffre 
et al. (2009). These simulations are conducted in a shearing box 
with dimensions (Jia, 4ffo, Hq'). As a check on the robustness of 
the main results and to verify their independence on the numer- 
ical implementation, some simulations were also performed with 
NIRVANA (Ziegler 1999, Papaloizou & Nelson 2003) adopting a 
shearing box with dimensions (L^, Ly, Lz) = {Hq, nHo, Ho)- 
The reference scale height is given by 



Ho 



CsO 



(UToV^^ 



(16) 



and thus corresponds to some prescribed isothermal sound speed 
CsO or, equivalently, temperature Tq. Note that as the simulations 
considered here are not isothermal, Cao need not necessarily cor- 
respond to the actual sound speed at any location or any time. In 
practice, however, it often approximates a volume average of the 
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initial average sound speed (for ZEUS runs) or the sound speed 
under quasi-steady conditions (NIRVANA runs). 

The box is periodic in the x and z direction and shearing pe- 
riodic in the azimuthal y direction (see Hawley et al. 1995). Space 
is scaled by Hq and time is measured in units of . Density is 
measured in terms of the average mass density in the box po, and 
magnetic field by the background field strength Bo (see below). 
Finally, temperature is scaled by To. 

The NIRVANA simulations employ a resolution of 
{N^,Ny,N^) = (128,200,128), where Ni is the number 
of grid cells in the i'th direction. This resolution level was adopted 
in Fromang et al. (2007) and Fromang (2010) and was there 
found to be converged. The ZEUS simulations, on the other hand, 
adopt a resolution of (128, 512, 128). Though more demanding, 
because of the higher resolution in the direction of shear, the 
ZEUS simulations can take advantage of the numerical benefits of 
an isotropic grid (Lesaffre & Balbus 2007). Additional runs were 
performed with the resolution reduced by a factor of two for both 
codes and these indicated no significant differences from the main 
results (see Section 5.4). 

In order to study simulations with a range of strengths of tur- 
bulent activity, three different magnetic field configurations (and 
boundary conditions) were implemented: fields with zero net-flux 
(or zero vertical-flux), fields with net vertical-flux, and fields with 
net toroidal-flux (see Hawley et al. 1995). In the net field cases, we 
must stipulate the initial plasma beta of a run: 



StvPq 
^0 



(17) 



where Bo is the strength of the mean field penetrating the box 
and Po — poHoti^ is the initial pressure. As /3 decreases, simu- 
lations become more active and display larger effective values of 
a. Furthermore, the cases with net vertical and with net toroidal 
flux produce turbulence manifesting similar mean values of a but 
with different characteristics (Hawley et al. 1995). This can give an 
indication of to what extent specification of the mean value of a is 
sufficient to characterize the thermal equilibria. 

In the zero flux case, we set the initial vertical field to be of 
the form Bq sm{x / Hq). However, provided Bo is large enough for 
the MRI to be resolved, the system attains a state after some 10 
orbits that is independent of the initial condition (Balbus & Hawley 
1998). 

For the diffusivities, we adopted parameters similar to those in 
Fromang et al. (2007) which compared results for zero flux simu- 
lations obtained from three codes, ZEUS, NIRVANA and the PEN- 
CIL code. These were found to be consistent and, in subsequent 
work of Fromang (2010), converged. We remark that Fromang et 
al. (2007) showed that the level of MRI turbulence critically de- 
pends on the magnetic Prandtl number Pm = v/rj, i.e. the ratio of 
kinematic viscosity to magnetic diffusivity. As a consequence, we 
set Pm = 4 to ensure sustained activity, and take 



= 32 X IQ'^Hln and -q ■- 



(18) 



Finally, in the case of NIRVANA, the initial data was prepared 
from similar isothermal simulations to those carried out in Fromang 
et al. (2007). Heating and cooling is then 'switched on' at a given 
time after a quasi-steady turbulent state is attained. The ZEUS runs 
begin from large amplitude perturbations which blend a large num- 
ber of MRI shearing waves, as calculated from the linear theory 
(Balbus Sl Hawley 1992). In this way a saturated turbulent state 



can be achieved relatively quickly and without an initial disruptive 
transient (Lesaffre et al. 2009). 



4.1 Numerical implementation of the radiative cooling model 

In this paper we are interested in thermal equilibria generated by 
models in which heating is provided directly through MRI turbu- 
lence. Accordingly, the cooling prescription must be implemented 
as an algorithm within our simulations. The values for Tc that are 
input into the function A (see equations l llOt and l |A13l l) is taken to 
be the vertically averaged temperature in the box. 



1 f^" 
Tc{x,y,t) = — T{. 

Jo 



x,y,z,t) dz, 



(19) 



and thus A is a function only of x, y, and time. Consequently, the 
cooling rate is the same for every z at a given horizontal location 
{x,y). This is a natural consequence of using a one-zone cooling 
model of the type we have adopted. 

As mentioned earlier, the governing equations that we numer- 
ically solve are scaled according to fiducial values for time, space, 
density, etc. However, when we evaluate the cooling function A, 
we must convert back to physical dimensions. This means we need 
to specify explicit values for Q, To, and E. All calculations were 
carried out adopting an angular velocity Q corresponding to a Ke- 
plerian circular orbit of radius 3 x 10^" cm around a star of mass 
IMq. This choice corresponds to the outer regions of an accre- 
tion disc around a white dwarf and yield Q — 2.22 x 10"'^ s~^ 
(FLP83). We explored, however, a range of values for the temper- 
ature scale To and E, over separate simulations, in order to sketch 
out the various thermal equilibria in the phase space of these vari- 
ables. The values were chosen so that the simulation begins near an 
estimated point on an S-curve calculated from an alpha-disc model. 
Once A is calculated in dimensional units it is rescaled to simula- 
tion units and fed into the energy equation In most cases, sim- 
ulations remain in the same phase space neighbourhood. However 
departures can occur in cases where there is a thermal runaway (i.e. 
no nearby equilibrium). In some cases the box dimensions become 
mismatched to the physical conditions, at which point the simula- 
tion was terminated. 

While To and E varied considerably, we employed two pa- 
rameter sets for {fi,E,X). The first set, corresponds to: fj, = 1, 
E = 3.5, A = 5. We refer to these as 'Set 1'. The second set was 
that adopted by FLP83: p = 0.5, E = 5.66, and A = 1. We refer 
to these as 'Set 2'. We choose two sets to extend the range of pos- 
sible disc states and also display the insensivity of our qualitative 
results to these parameters. 



4.2 Diagnostics 

Of central importance in our simulations are the turbulent stresses 
which, via the transport of angular momentum, facilitate mass ac- 
cretion and localised heating (see Balbus & Papaloizou 1999). We 
measure the stresses via the turbulent alpha: 



{Txy) _ {pVxVy — BxBy/{4:n)) 

(P) ~ (P) 



(20) 



where Txy is the xy component of the total turbulent stress ten- 
sor and the angle brackets indicate a box average. This is a quan- 
tity that fluctuates in time, and thus must be distinguished from 
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Figure 2. A screen shot of the temperature in Kelvin during a saturated 
equilibrium state. Pai'ameter Set 1 is used with S = 1000 glcw? with net 
toroidal field. The upper branch solution has been achieved. 



the (constant) Shakura-Sunyaev alpha parameter ass that is used in 
laminar disc modelling, defined through 



(21) 



The two quantities should coincide once the turbulent alpha is av- 
eraged over long times. 

The turbulent stresses Txy extract energy from the orbital 
shear and, once this energy travels down a turbulent cascade, it is 
thermalised by Ohmic and viscous dissipation. This turbulent heat- 
ing is captured directly in our model via the first two terms on the 
right side of (|6}. In a thermal quasi-steady state this heating rate eo 
must balance the radiative cooling rate A. If this is not achieved the 
temperature of the system T will evolve until a thermal balance is 
obtained. 

In a steady state laminar disc the accretion rate M is given by 



(22) 



where v is the vertically averaged value of v weighted by density. 
In a turbulent disc, the equivalent relation employing local averaged 
quantities is 

M = y {Tr^)dz ^ 2-kHq i-^, (23) 
and, like a, is a time varying quantity. 



5 NUMERICAL RESULTS 

We undertake numerical simulations of MRI-induced turbulence 
for two sets of cooling parameters. For each parameter set we in- 
vestigate three different magnetic configurations: (a) zero-net flux, 
(b) net-vertical flux, and (c) net-toroidal flux. And for each con- 
figuration we conduct a suite of separate simulations with various 
(S, To) in order to sketch S-curves in each of the 6 cases. 

In the net-vertical flux runs, the initial plasma beta is /3 = 10* . 
This large value is taken in order to minimise recurrent channel 
flows — an artefact of small boxes — that may skew our equilib- 
rium results (Hawley et al. 1995, Sano & Inutsuka 2001, Latter et 




100 200 300 400 500 600 

Surface Density (g/cm^) 

Figure 4. This figure summarises the simulations of Fig. |3]by plotting ei- 
ther their time averaged M in the (E, M) plane, if a thermal equilibrium 
is achieved, or by showing the initial and end M if the simulation exhibits 
a thermal runaway. Equilibria are plotted as dots with the eiTor bars indi- 
cating the range of fluctuations. Runaways are plotted as arrows. Note that 
despite the range of the fluctuations, equilibria (dots) never flipped between 
branches. Only simulations begun near branch edges transitioned. The an- 
alytic S-curve (dashed), provided for reference, is obtained from following 
the procedure of FLP83 for ass = 0.034. 



al. 2009). For the net-toroidal flux, stronger fields were employed 
and the initial beta was either /3 = 50 or 100. 

Once heating and cooling were activated, each simulation was 
run until the system had relaxed towards a turbulent thermal equi- 
librium, if one was nearby, or had monotonically heated up or 
cooled down so as to change solution branches (cf. the vertical ar- 
rows in Fig. 1). In the former cases, the relaxation time depended on 
the heating rate and hence the magnitude of the turbulent stresses. 
Because of the small values of ct in the zero flux runs, this process 
could take some time and these runs were typically continued for 
between 150 and 250 orbits. The net toroidal-flux runs exhibited 
larger a and a corresponding shorter heating time scale. They re- 
quired 30 - 100 orbits to reach a quasi-steady states with meaning- 
ful statistics. Fig.|2]shows a screen shot of the turbulent vertically- 
averaged temperature of a representative run once a thermal equi- 
librium has been achieved. Here a net-toroidal flux simulation has 
been selected and the upper branch achieved. 

In the next subsection we present a collection of 8 runs asso- 
ciated with parameter set 1 and a net-toroidal flux. We show their 
thermal solution trajectories as functions of time and how these tra- 
jectories change as we vary the surface density E. In this way we 
can trace out individual thermal states along a limit cycle in the 
(E, M) plane. Our full set of numerical results are then displayed 
with a more detailed discussion in the following subsection, where 
S-curves are sketched out for all magnetic configurations and pa- 
rameter sets. These, however, are presented in terms of E and box- 
averaged temperature. Once this is done we discuss the weak de- 
pendence between the turbulent stresses and the heating, and finally 
issues of numerical convergence. 
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Figure 3. The time evolution of the mass accretion rate M in eight simulations for the case of net toroidal-flux and parameter set 1. Here /3 = 100. Each 
simulation is distinguished by its (conserved) surface density S in cgs units. In a sequence running from left to right and top to bottom, the conserved surface 
density moves first through increasing values for which there is a mean steady state corresponding to the lower branch (LB) of the S-curve plotted in Fig.O 
The simulation that is second from the top and on the right then moves to the upper branch (UB) as there is no available mean steady state on the lower 
branch. The sequence resumes with states of decreasing surface density tracing out the hot high-accreting solution on the upper branch. The final bottom right 
simulation then transitions to the lower branch, there being no steady state available on the upper branch for its low E. 



5.1 Simulation tracks in the (E, M) plane 

In laminar disc modelling it is common to describe local S-curves 
and limit cycles within the phase plane of E and M because it is the 
accretion rate M that dominates the state of the disc. Accordingly, 
we begin by discussing a subset of our simulation results in terms 
of the evolution of M. 

In Fig. [5] we plot the evolution tracks of eight simulations with 
net-toroidal flux, with j3 = 100, and parameter set 1. The sim- 
ulations differ in their (conserved) surface density E and initial 
temperature To. The simulations achieve either a thermal quasi- 
steady equilibrium near their initial state, in which M fluctuates 
about a well-defined mean value, or they catastrophically heat up 
or cool down from their initial state, with an accompanying in- 
crease/decrease in M. These results are summarised by Fig.|4]in 
the (E, (A/)) plane, where the (A/) of each simulation is com- 



puted from a time average over the course of the run. Runs that 
achieved a thermal steady state near their initial condition are rep- 
resented as dots with error bars showing the range of the fluctua- 
tions in M. Runs that evolved significantly away from their initial 
condition are represented by arrows, the base of which denote their 
starting point and the tip of which their end point when the simula- 
tion was terminated. 

As we progress from left to right and top to bottom through 
the eight panels in Fig. [3] we track the S-curve in Fig. |4] starting 
at the bottom left comer. First we travel towards the right along the 
lower branch, then up to the higher branch, then to the left along the 
upper branch, and finally back down to the lower branch. For the 
first three simulations/panels, the surface density E moves through 
increasing values, each of which permits a cool low-accreting equi- 
librium. But the fourth simulation rapidly migrates from the vicin- 
ity of the lower branch: no cool low-accreting steady state is avail- 
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able, because it possesses a E that is too large. It instead heats up 
and approaches an upper branch of hot solutions characterised by 
efficient accretion. The sequence then resumes with a procession of 
hot equilibria with decreasing E. The eighth simulation transitions 
from the vicinity of the upper branch to the lower branch as there 
is no hot high-accreting state available that corrresponds to its low 
surface density. 

Though a single run in local geometry cannot describe a com- 
plete limit cycle from beginning to end, our sequence of multiple 
runs can describe the constituent parts of a cycle, state by state. 
This is permitted because of the separation of time-scales between 
the fast turbulent/thermo dynamics that determine the equilibria 
(> 1/0, ~ 10^ s) and the slow dynamics of the cycle, which is 
governed by the accretion rate (> l/{aO) ~ 10^ s). 

Despite the large range of the M fluctuations witnessed in 
Fig.|4]there is never any danger that the equilibrium states sponta- 
neously undergo transitions. States on the lower and upper branches 
are robust and distinct, with their stability tied to (the smaller) vari- 
ations in T rather than M. This point is is more apparent when the 
S-curves are plotted in the (E, T) plane (see next subsection). 

Finally, we have superimposed on Fig.|4]an analytic S-curve 
computed using the formalism of FLP83 with a Shakura-Sunyaev 
alpha of Qfss = 0.034. The location of the two numerical branches 
and the stable analytic branches is reasonably consistent. However, 
the numerical branches often seem more extended to the right. 

5.2 Thermal equilbria: S-curves as viewed in the (E, (T)) 
plane 

In this subsection we present the bulk of our simulation results 
and plot S-curves for the cases of both parameter sets and for all 
three magnetic field configurations. The same qualitative behaviour 
emerges in every scenario, which reinforce the robustness of these 
features. 

Tables [T] and |2] summarise the simulation results and set-up 
for the runs that quickly established thermal equilibria. Simulations 
which underwent a thermal evolution, corresponding to either sec- 
ular heating or cooling, are described in Tables [3] and |4] In Figs|5] 
and|6]we plot S-curves for the various runs. These show the time- 
averaged temperature (once thermal equilibrium is achieved) as a 
function of surface density E (a conserved input parameter). The 
figures graphically summarise some of the data in Tables 11141 As 
in Fig. 4, the blue points with error bars represent the runs which 
were begun near a stable thermal equilibrium. The arrows repre- 
sent 'thermally unstable' runs, which were begun at the base of 
each arrow and were evolved to a state at the tip of each arrow, 
at which point the simulation either approached a stable branch or 
was stopped. 

First, the simulations reveal that the volume averaged temper- 
ature fluctuations in the turbulent equilibria are relatively small, 
some 10%, much smaller than for M. This means that the ther- 
mal equilibria are in fact very stable: turbulent fluctuations are not 
so extreme as to push the system over the unstable intermediate 
branch in the (E, T) diagram. In all the cases we simulate, equilib- 
rium states never jump branches spontaneously. Only at the 'cor- 
ners' of the S-curve is this a prospect. And on each corner the sys- 
tem could only move in one direction. As a consequence, limit cy- 
cle behaviour should follow broadly along the lines predicted by 
the classical theory. 

The 'runaway' simulations represent the system switching 



branches on account of the fact that no equilibrium state was avail- 
able. The simulations that went from the vicinity of the lower to 
upper branches heated up at a rate afi, as expected. The simula- 
tions that left the vicinity of the upper branch cooled at a varying 
rate determined by A. In some cases turbulence in the latter simu- 
lations died out once the volume averaged temperature in the box 
achieved a level much lower than the starting temperature Tq. This 
is because the scale height in the box becomes so small (of order the 
dissipation scale) that the MRI is suppressed. Similarly, the physi- 
cal condition of the gas becomes mismatched to the box size when 
the simulations catastrophically heat up. Simulations in both cases 
were discontinued once this occurred. 

Local simulations that begun near the S-curve 'corners' would 
often bobble around indeterminately before switching curves. In 
these cases the system was near criticality and therefore small fluc- 
tuations in the heating rate change the stability of the system un- 
ceasingly. A sudden change in heating may cause the local equi- 
librium to vanish; the gas heats up or cools only for the heating 
to readjust and the equilibrium state to reappear. It follows that in a 
limit cycle the branch jumping near branch ends can have a stochas- 
tic component. However, we reiterate we do not see such transitions 
between the central regions of the upper and lower branches. 

Superimposed upon the numerical S-curves are representa- 
tive analytical S-curves derived from FLP83, in which ass is a 
user defined parameter. These ass values were chosen so that good 
agreement was obtained with the numerical equilibria on the lower 
branch. Note, in particular, that once we set the analytic curves to 
match the lower equilibria, the upper numerical equilibria tend to 
exhibit slightly higher temperatures than expected from the analytic 
S-curves. Quite generally, the upper numerical branch is better ap- 
proximated by alpha models with larger ass values than required 
on the lower branches. This could reflect the fact that the effect 
of turbulent heating is not accurately described by the classical a 
model given the run times of the simulations. However, it could 
also to some extent reflect a bias in the way the simulations were 
set up and carried out. 

A close study indicates that the turbulent a weakly depends 
on the ratio of the average temperature (T) to the reference tem- 
perature To, or equivalently the mean scale height to the reference 
scale H(i . Very 'hot' (small) boxes with T > To exhibit a smaller a, 
while 'colder' (large) boxes yield larger a, with relative differences 
~ 30%. This is essentially a numerical effect arising from varia- 
tions in the effective Reynolds number and it can operate within 
any simulation. As a consequence, there is some low level indeter- 
minacy in both the properties of the equilibrium solutions and the 
exact locations of the S-curve corners. Their existence and gross 
properties, however, remain. Another, less important, consequence 
of this 'alpha feedback' is the larger amplitude and intermittent 
thermal fluctuations in the saturated turbulent states, as compared 
to isothermal models. 

We remark that the need for a branch-dependent a, in the con- 
text of a disc modelling, was recognised early in order to obtain 
decent matching of outburst behaviour with observations (Smak 
1984a, Lasota 2001). However, for the reasons mentioned above, 
caution should be exercised in the interpretation of our numerical 
local a determinations in this context. Finally, note that the levels 
of turbulence and values of a increase with imposed toroidal or 
vertical magnetic flux. Thus if the flux were to increase as a system 
transitioned from a low state to a high state, the higher state would 
be associated with a larger value of a. Thus global disc simulations. 
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Figure 5. S-curves computed with ZEUS for parameter set 1: (/i, E, A) = (1, 3.5, 5) and for the thi'ee magnetic configurations. The first panel coiresponds 
to zero net-flux, the second panel to net vertical-flux (initial /3 = 10"*), and the third panel to a net toroidal-flux (initial /3 = 100). Simulation results are 
indicated on the (S, T^) plane, with T^, in K, evaluated as the volume average of the temperature, T, over the box. E is expressed in cgs units. Runs that 
attained an approximate quasi-steady state are indicated by a dot, corresponding to the time average mean volume averaged T, and with errorbars, indicating 
the maximum and minimum volume averaged T once the steady state is achieved. Cases indicated with arrows pointing upward/downward underwent a 
heating/cooling instability, with the volume averaged T on average showing a monotonic increase/decrease with time. Notation is otherwise the same as 
in Fig. O The dashed curves represent equilibrium S-curves calculated using the formalism of FLP83 for ass = 0.007, 0.0275, 0.034 for the zero flux, 
net-vertical flux, and toroidal flux cases respectively. 



for which net-flux need not be conserved locally, could potentiaOy 
lead to such local a variations. 



5.3 Turbulent stresses and heating 

We compared the box averaged pressure (and temperature) with 
the turbulent stress in our quasi-steady equilibria. In Fig. [7] we plot 
a representative evolution of these two fluctuating quantities. As is 
clear, on the orbital time-scale there only exists a weak dependence 
between pressure and the turbulent stress (described by a). There 



are three discrepancies. First, pressure lags behind the turbulent 
stress by roughly an orbit, which indicates that on short times the 
causal relationship is from viscous stress to pressure and not the 
other way around, in agreement with Hirose et al. (2009). Peaks 
in a describe intense large-scale events, the energy of which sub- 
sequently tumbles down a cascade to reach dissipative scales. The 
gas heats up and the pressure increases. This process, however, is 
not instantaneous and the original a signal is significantly modi- 
fied ('smoothed-out') by the time that it manifests itself as heating 
(see Pearson et al. 2004 for similar results in forced isotropic tur- 
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Figure 6. S-curves computed with NIRVANA and ZEUS for the second parameter set: E, A) = (0.5, 5.66, 1). The first panel corresponds to zero net- 
flux, the second panel to net vertical-flux (initial /3 = 10*), and the third panel to a net toroidal-flux (initial j3 = 50). The dashed curves represent equilibrium 
S-curves calculated using the formalism of FLP83, with ags = 0.009, 0.0475, 0.07 for the zero flux, net-vertical flux, and toroidal flux cases respectively. 
Notation is otherwise the same as in Fig. [3] 



bulence). Second, the volume averaged pressure and a time series 
are dissimilar; in particular, the volume averaged pressure appears 
at best as a heavily convolved version of a. Lastly, variations in 
a are much larger than related variations in the volume averaged 
pressure, which are relatively mild, with the former ~30%, and the 
latter less than 5%. 

This poor correlation between the viscous stress and pressure 
on short times casts doubt on the existence of thermal instabili- 
ties that rely on these two quantities to be functionally dependent 
(Shakura & Sunyaev 1976, Abramowicz et al. 1988). Indeed, nu- 
merical simulations of the MRI in radiation pressure dominated 
discs reveal no thermal instability (Hirose et al. 2009). Recently, 
linear stability analyses have been conducted that allow the pres- 
sure to lag behind stress, and vice-versa. These show that, de- 



spite this impediment, instability can in fact persist, especially with 
the short time-lags witnessed in our simulations (Lin et al. 2011; 
Ciesielski et al. 2012). It would appear then that the main obsta- 
cle to instability in turbulent models of discs is due to something 
else, perhaps the stochastic nature of the viscous stress which is 
not faithfully reproduced in the behaviour of the volume averaged 
pressure (Janiuk and Misra 2012). 



Finally, note that over long times the turbulent stress Txy and 
pressure are correlated. Thus turbulent accretion on the hot upper 
branch of equilibria is greater than accretion in the colder states (cf. 
Section 5.1). 
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Table 1. Summary of the inputs and average values of simulations undertaken by ZEUS with parameter set I: fj. = 1, E = 3.5, A = 5. The magnetic field 
configuration is varied as is the input surface density S. The initial f3 is 10* for the net vertical field runs, and 100 for the net toroidal field runs. Quasi-steady 
thermal equihbria are attained that align with one of the two stable branches of the S-curve, with 'l.b.' denoting 'lower branch', and 'u.b.' denoting 'upper 
branch'. 



5.4 Numerical convergence 

Finally, we briefly discus.s is.sues associated with the numerics. Pa- 
rameter set 1 was primarily undertaken with ZEUS and parameter 
set 2 with NIRVANA, however the two codes were run on an over- 
lap set of some 8 simulations. The comparison yielded good agree- 
ment in both (T) and a, with discrepancies within the range of the 
fluctuations. It follows that not only are our qualitative results ro- 
bust with respect to parameter choices and magnetic configurations, 
they are robust with respect to numerical method. 

In addition, we performed a short study in order to confirm 
that the simulations were adequately converged with respect to 
resolution and Reynolds nvunber. A representative sample of runs 
were taken and we evaluated the quantities a and (T) for different 
Reynolds numbers while keeping the Prandtl number constant and 
equal to 4. In summary we found that doubling the kinematic vis- 
cosity f (halving the Reynolds number) led to changes in a and (T) 
that were less than 5%. This reproduces the very weak Reynolds 
number dependence uncovered by Fromang (2010). Moreover, the 
result held whether we took 64 grid points per Ho or 128 points. 



Thus we were assured that these quantities were satisfactorily con- 
verged. 

Next we kept the diffusion coefficients fixed and varied the 
resolution. We evaluated the above quantities at both 64 and 128 
grid cells per Ho. The lower resolution runs yielded slightly lower 
(T) but the discrepancy was within 10%. This small decline is un- 
derstandable because at coarser resolutions more of the turbulent 
energy is removed by the grid and thus less captured by physi- 
cal Ohmic and viscous heating. As a consequence, there is slightly 
less direct heating. That said, the effect is small and in general no 
greater than the statistical fluctuations of the runs. It leads in no 
way to any change in the qualitative behaviour exhibited. 



6 CONCLUSION 

In this paper we performed a suite of numerical simulations of the 
MRI in local geometry with the ZEUS and NIR"VANA codes. Both 

viscous and Ohmic heating is included, while the radiative cool- 
ing is approximated by a physically motivated cooling function 
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4.67 X 10" 


3.31 X 10" 


0.041 


40 


Net Vertical (u.b.) 


200 


4.33 X 10 


A T c >. * 1 r\4 

4.75 X 10 


3.9U X 10 


0.028 


59 




250 


5.7 X 10'' 


6.12 X 10 


5.28 X 10 


0.038 


25 




3UU 


D.U X iO 


6.34 X 10 


5.66 X 10 


0.042 






80 


4.20 X 10^ 


4.,33 X lO-'' 


4.07 X 10^ 


0.061 


38 


Net Toroidal (l.b.) 


100 


4.31 X 10^ 


4.43 X 10-^ 


4.20 X 10^ 


0.068 


50 




130 


4.41 X 10^ 


4.54 X 10-^ 


4.30 X 10^ 


0.071 


34 




170 


4.89 X 103 


5.03 X 103 


4.68 X 10^ 


0.08 


44 




200 


5.47 X 10^ 


5.87 X 10^ 


5.10 X 10^ 


0.072 


44 




100 


3.31 X 10* 


3.63 X 10" 


3.02 X 10" 


0.085 


43 




140 


4.49 X 10" 


4.69 X 10" 


4.24 X 10" 


0.075 


45 


Net Toroidal (u.b.) 


160 


4.84 X 10" 


5.20 X 10" 


4.58 X 10" 


0.072 


37 




200 


5.18 X 10" 


5.73 X 10" 


4.88 X 10" 


0.054 


37 



Table 2. Summary of the inputs and average values for simulations undertaken by NIRVANA adopting parameter set 2: fi = 1, E = 5.66, A = 1. The initial 
/3 is lO" for the net vertical field runs, and 50 for the net toroidal field runs. Some simulations are performed with the same value of S but different initial 
mean temperatures: * indicates a run begun wifii (T) = 3.85 X lO", whereas ** indicates (T) = 4.75 x lO"; f denotes an initial (T) = 5.20 x lO" and ft 
denotes (T) = 6.00 x lO". The larger variations associated with the former pair of cases occurs because the steady state thermal equilibrium is expected to 
be located near the upper bend in the S curve. 



that summarises the strong effect of temperature on the ability of 
the disc gas to retain heat. Different magnetic configurations (zero 
flux, net-toroidal flux, net-vertical flux) and parameters were tri- 
aled, with little change in the qualitative results. 

Our simulations unambiguously exhibit the development of 
thermal instability and hysteresis. In particular, through a sequence 
of runs we can sketch out characteristic S-curves in the phase space 
of (S, Tc) and (E, M), which are central to the classical outburst 
model. It hence appears that MRI turbulence is not so intermittent 
as to endanger the robustness of the cycle. Temperature fluctuations 
are well-behaved and relatively small and there is no spontaneous 
jumping from one stable branch to the other. Only near the 'cor- 
ners' of the S-curve does significant stochasticity enter, as then the 
existence or not of a local equilibrium is uncertain. This feature 
will add some degree of low level 'noise' to the observed outburst 
time-series. In addition, the a we record on the two stable branches 
differ slightly but systematically. Because of the constraints of our 
local model, this result is more suggestive than anything else. How- 



ever, it does indicate that in global disc simulations we may indeed 
see a systematic difference in the two branches, as required by the 
classical theory. 

Finally, on the orbital time the turbulent stresses and pressure 
only weakly depend on each other in our simulations. Pressure al- 
ways lags behind the turbulent stress, and thus the causality is from 
the stress to the pressure, via the turbulent heating (in agreement 
with Hirose et al. 2009). Moreover, the pressure response is a sig- 
nificantly 'smoothed out' echo of the turbulent stress features. Both 
facts indicate that thermal instability driven by turbulent heating 
variations, in which the stress is a function of pressure, may not 
operate in real discs. On long times, however, there is necessarily 
a feedback of the pressure on the stress, which leads to different 
accretion rates on the two branches of the S-curve. 

Our work presents a first step towards unifying simulations of 
full MHD turbulence with the correct thermal and radiative physics 
of outbursting DNe and LMXBs, and possibly young stellar ob- 
jects. We have begun with with the most basic model that yields 
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Parameter Set 1 : Thermal 


runaways 




B Configuration 


E 


initial T 


final T 


Orbits 




300 


3.7 X lO* 


^^ n v 1 n3 
O.U X lU 


ou 


Zero Vertical (c) 


400 


4 X 10'^ 


o.O X lU 


1 




500 


4.75 X lO" 


2.7 X 10^ 


178 


Zero Vertical (h) 


1700 


4.5 X 10^ 


3.4 X lO'' 


535 




150 


4.8 X lO"* 


6.3 X 10^ 


20 


Net Vertical (c) 


175 


6.5 X lO** 


7.1 X 10^ 


31 




188 


6.0 X lO"* 


7.8 X 10^ 


52 


Net Vertical (h) 


600 


3.3 X 103 


1.1 X lO'' 


101 


650 


3.3 X 10« 


1.4 X lO"* 


108 


Net Toroidal (c) 


100 


5.0 X lO"* 


9.0 X 10^ 


14 




560 


3.6 X 10^ 


2.5 X 10* 


43 


Net Toroidal (h) 


625 


3.6 X 10^ 


4.1 X 10'' 


73 




750 


3.0 X 10^ 


4.4 X 10"' 


68 



Table 3. Summary of the inputs and range of values found for simulations 
undertaken ZEUS adopting parameter set 1: /.t = 1, E = 3.5, A = 5. 
The magnetic field configuration is varied as is the input surface density S. 
These simulations resulted in thermal runaways, i.e. persistent heating or 
cooling as indicated by either (h) or (c) . This is because they were initially 
located near the corners of S curves in the surface density time averaged 
mean temperature plane. In the latter case this occured on a time scale of a 
few orbits in some simulations. In general very large changes to the mean 
temperature in these simulations eventually resulted in the computational 
box becoming mismatched to the putative scale height. Accordingly, they 
were then not continued. 



Parameter Set 2: Thermal runaways 



B Configuration 


E 


initial T 


final T 


Orbits 


Zero Vertical (c) 


100 
186 


3.18 X 10* 
3.18 X 10* 


3.18 X 10^ 

3.19 X 103 


5.6 
17 




686 


8.49 X 103 


1.70 X 10* 


270 




850 


6.34 X 10^ 


1.78 X 10* 


155 




950 


4.75 X 103 


1.50 X 10* 


54 


Zero Vertical (h) 


1000 


6.37 X 10^ 


4.07 X 10* 


119 




1100 


2.00 X 10^ 


1.17 X 10* 


92 




1200 


6.00 X 10^ 


3.18 X 10* 


254 




1200 


8.00 X 10^ 


5.92 X 10* 


215 


Net Vertical (c) 


80 
100 


3.82 X 10* 
3.18 X 10* 


6.37 X 10« 
5.31 X 10^ 


3.8 
5.25 


Net Vertical (h) 


275 
300 


2.00 X 10^ 
1.27 X 10* 


1.62 X 10* 
1.78 X 10* 


81 
33 


Net Toroidal (c) 


80 


3.82 X 10* 


1.59 X 10* 


5.9 


Net Toroidal (h) 


225 
250 


4.00 X 10^ 
9.55 X 10^ 


5.56 X 10* 
1.46 X 10* 


66 
43 



Table 4. As in table [3] but for simulations undertaken by NIRVANA and 
ZEUS adopting parameter set 2: fj. = 0.5, E = 5.66, A = 1. 




Time (orbits) 



Figure 7. Plots of the turbulent stress, quantified by a, versus time and 
the volume averaged pressure versus time. Parameters correspond to Set 1 ; 
there is a net toroidal-flux associated with initial /? = 100, and S = 250 
(UB) (see Table 1). It is clear that the latter tracks the sfi'ess imperfectly, 
with much of the short time scale stochastic variabiUty seen in the stress 
washed out. Gross features are also associated with a time-lag of about an 
orbit. 

the correct physics: local simulations of gas inhabiting a single 
radius in the disc with a simple radiative prescription. This is a 
'ground zero' test of the compatibility of the MRI with the puta- 
tive thermal limit cycles of outbursting discs. If MRI-turbulence 
had failed in this basic setting it would probably fail in more ad- 
vanced models as well. Now that we are assured of this compatabil- 
ity, a variety of further work may be attempted. For instance, sim- 
ulations could be performed in vertically stratified shearing boxes 
with more realistic radiation physics (in the flux-limited diffusion 
approximation with appropriate opacities). In addition, cylindrical 
MRI simulations should be performed with the FLP83 cooling pre- 
scription, anologous to the global a disc calculation of Papaloizou 
et al. (1983). Such simulations would describe how real turbulence 
mediates the heating and cooling fronts that propagate through the 
disc during a transition between branches. 
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APPENDIX A: RADIATIVE COOLING MODEL 

In this appendix we detail how the coohng prescription of Section 
3.2 was derived. In particular we show how the effective tempera- 
ture Te is calculated within each of the three important gas regimes 
introduced. 



Al Regime 1: the optically thick hot regime 

The radiative cooling function can be written as a flux in the diffu- 
sion approximation by A = V ■ f where 



3Kp 



VT, 



(Al) 



in which k is the opacity, a — 4cr/c is the Stefan-Boltzmann radi- 
ation constant and c is the speed of light. Setting the optical depth 
as 



i^{x,y,z')p{x,y,z')dz' , 

we obtain the vertical component of F in the form 

p 4acT^ dT 



3 di 



(A2) 



(A3) 



We integrate from near the disc surface where T — Te, (the effec- 
tive temperature), and r = Te ~ 1 to the mid plane where 2 = 0, 
T — Tc and t — Tc- Thus we obtain 



(A4) 



where Tc is the mid plane temperature. 

In the hot optically thick Regime 1, ^ Te ~ 1 and T^ S> 
Te so that l lA4b implies 



T^ 



(A5) 



where we have approximated to be constant and equal to its sur- 
face value acTc /A. This equation relates the effective temperature 
Te to the central temperature Tc- 

To complete the prescription we take the scale height to be 



(A6) 



and set the surface density to 

E = 2pcH. 



(Al) 



Finally, as most of the optical depth arises from the regions near the 
midplane, we approximate Eq. l lA2t by 



Tc — KcE/2, 



(A8) 



where Kc is the opacity evaluated at the density pc and the midplane 
temperature Tc. In this hot ionised regime, we can approximate k 
by the following formula: 



K = 1.5 X lo^v^"^ ■^ 



(A9) 



(FLP83). So for specified fi, the above relationships enable Tc (and 
hence A) to be related to conditions in the mid plane, and hence Tc 
and E. 



A2 Regime 2: the warm transitional regime 

This regime exhibits a cooler disc midplane that is still well ionised 
but surface layers that are < 5000 K, and which are poorly ionised. 
As a result, near the photosphere k not only drops significantly but 
becomes a steeply increasing function of temperature. It follows 
that the requirement t ~ 1 at the photosphere determines the ther- 
mal structure of the disc, as in cool stars (Hayashi & Hoshi 1961). 
Instead of iA5t one must impose the condition Te ~ 1, where Te 
is the optical depth at the location where T — Tc. As the optical 
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surface is ~ above the mid plane, it converts to a relationship 
between the disc surface density Te and Tc that leads to the central 
unstable part of the S-curve when thermal equilibria are considered. 

We introduce the quantity r* which we define through r* = 
Ke S with Ke being the opacity evaluated at the central density pc 
but with the photospheric temperature T^. It is anticipated that r* 
overestimates the optical depth of material above the disc photo- 
sphere Te by some factor of order unity. This factor we quantify 
through the dimensionless constant E, so that 

Erl^Te-l (AlO) 

Once the dependence of on Te and pc is specified, equation 
l lAlOl l supplies a means to obtain T^, and hence A, when the disc 
is in Regime 2. Finally, an approximate functional form of k in this 
regime of incomplete hydrogen ionisation is 

K = lG-'"'p^''^T^\ (All) 

(FLP83). 

A3 Regime 3: the cool optically tliin regime 

In this regime the vertical structure of the disc is isothermal with 
T ~ Tc. However, because the disc is optically thin in the vertical 
direction, Fz is reduced from acT* /4 by a factor ~ 2rc. Thus we 
set 

= 2\t,T^, (A12) 

where A is a dimensionless constant of order unity that takes ac- 
count of appropriate frequency averaging and any other necessary 
corrections arising from a more complete discussion of radiation 
transport. This condition also connects Te to Tc and E. The central 
optical depth in this regime Tc is computed using Eq. 1 IA8I 1 and the 
opacity prescription introduced in Eq. JAIU . 



9000 




Figure Al. The effective temperature Te as a function of central temper- 
ature Te as determined from the FLP83 cooling prescription, cf. Eq. 
Parameters are p = 1, E = 3.5 and A = 5. In addition, we have indicated 
on the curve which regime is dominant. 



A4 Interpolated form 

The prescription for Te is given in Eq. jilt for each of the three 
regimes. In numerical simulations one could switch from one form 
to the other according to the local value of Tc; we however em- 
ploy an analytic interpolation of the three forms, as in FLP83. In 
Eq. Jilt , the effective temperature Te takes one of three forms, 
which we denote by either Ti, T2 or T3 where the subscript in- 
dicates the associated regime. These three functions are then inter- 
polated via the following: 

_ iTt^Ti)Ti ^^^3^ 



m+Tiy/^+Tif' 

(FLP83). It is this algebraic expression for Te that is used in Equa- 
tion l llOt to evaluate A. The two free parameters A and E can be 
chosen to fit results obtained from vertical integrations with stan- 
dard local Q modelling. In Fig. (Al) we plot the resulting func- 
tional dependence of Te on the central temperature Tc with the 
three regimes indicated on the curve. In regimes 1 and 3 the effec- 
tive temperature increases steeply with central temperature, while 
in the transitional regime 2 the dependence is very weak on account 
of the rapid changes in opacity near the surface. In this intermedi- 
ate regime the radiative cooling rate is very 'flat' and, consequently, 
thermal instability ensues (cf. Eq. 
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